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Abstract.  An  immersed  interface  method  for  solving  linear  elasticity  problems  with 
two  phases  separated  by  an  interface  has  been  developed  in  this  paper.  For  the 
problem  of  interest,  the  underlying  elasticity  modulus  is  a  constant  in  each  phase 
but  vary  from  phase  to  phase.  The  basic  goal  here  is  to  design  an  efficient  numerical 
method  using  a  fixed  Cartesian  grid.  The  application  of  such  a  method  to  problems 
with  moving  interfaces  driving  by  stresses  has  a  great  advantage:  no  re-meshing  is 
needed.  A  local  optimization  strategy  is  employed  to  determine  the  finite  difference 
equations  at  grid  points  near  or  on  the  interface.  The  bi-conjugate  gradient  method 
and  the  GMRES  with  preconditioning  are  both  implemented  to  solve  the  resulting 
linear  systems  of  equations  and  compared.  Numerical  results  are  presented  to  show 
that  the  method  is  second-order  accurate. 

Keywords:  elasticity,  interfaces,  jump  conditions,  finite  differences,  the  immersed 
interface  method. 
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1  Introduction 

Elasticity  problems  of  multiple  phase  elastic  materials  separated  by  phase  interfaces 
often  arise  in  materials  science  [18,  21].  Two  important  examples  of  such  problems 
occur  in  the  microstructural  evolution  of  precipitates  in  an  elastic  matrix  due  to  the 
diffusion  of  concentration  and  in  the  morphological  instability  due  to  stress-driven 
surface  diffusion  in  solid  thin  films,  cf.  e.g.,  [2,  10,  13]  and  the  references  therein. 
The  understanding  of  these  physical  processes  is  crucial  to  improve  material  stability 
properties,  and  in  turn  to  develop  new  and  advanced  materials  that  have  applications  in 
automobile  manufacture,  aircraft  industries,  and  modern  communication  technologies. 
However,  solving  such  elasticity  problems  are  often  very  difficult  due  to  complicated 
geometries,  multiple  components,  and  nonlinearities  that  appear  in  these  problems.  For 
these  reasons,  there  has  been  a  great  interest  recently,  in  all  materials  science,  scientific 
computing,  and  applied  mathematics  communities,  in  developing  efficient  and  accurate 
numerical  techniques  for  elasticity  problems  with  interfaces  separating  multiple  phases. 
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In  this  paper,  we  consider  a  two-dimensional  problem  that  can  arise  from  many 
modeling  situations  such  as  two-phase  elastic  plates  in  the  setting  of  plane  strain  or 
plane  stress  [5,  9,  20].  We  assume  that  the  elastic  material  occupies  a  bounded  domain 
Fl  C  K2  with  boundary  To  =  dFl.  For  our  computational  purpose,  we  assume  that  the 
domain  Fl  is  rectangular.  The  two  phases  of  the  material  occupy  regions  Fl+  and  fl-, 
respectively,  so  fl  =  fl+  U  fl-  and  fl+  fl  fl-  =  0.  We  denote  by  F  =  fl+  fl  fl-  the 
interface  separating  these  two  phases,  cf.  Figure  1.1. 


Figure  1.1:  The  geometry  of  the  underlying  problem. 


The  equilibrium  equation,  interface  conditions  on  T,  and  the  boundary  conditions 
on  Tq  [5,  9,  18,  20]  are: 


V  •  a  +  F  =  0 

in  Fl+  U  Fl~, 

(1.1) 

M  =  o 

on  r, 

(1.2) 

[an]  =  T 

on  r, 

(1.3) 

u  =  u0 

on  r0, 

(1.4) 

where  a  is  the  stress  tensor,  F  =  (Fi,  F2)1  :  fl  — >  M2  is  the  body  force  which  is  known 
(a  superscript  T  denotes  the  transpose),  u  =  {u\,U2)T  :  fl  — >  M2  with  u\  =  u\(x,y )  and 
U2  =  U2{x,y)  is  the  displacement  field,  for  a  function  v,  [V]  =  v+  —  v~  with  v±  =  v|q± 
denotes  the  jump  of  v  across  the  interface,  n  =  (n\.  712) 1  is  the  unit  normal  vector 
to  the  interface  T,  pointing  from  the  —  phase  to  the  +  phase,  T  =  is  a  given 

vector-valued  function  on  the  interface  F  which  measures  the  jump  of  the  traction 
cm  across  the  interface  T,  and  uo  is  a  given,  vector-valued  function  that  represents  the 
displacement  on  the  boundary  To-  The  jump  condition  (1.2)  means  that  u  is  continuous 
across  the  interface.  It  indicates  that  the  underlying  material  has  no  fracture.  Notice 
that  we  have  introduced  a  generally  nonzero  data  T  for  the  jump  of  the  traction  in 
(1.3).  This  will  be  useful  to  model  more  general  physical  situations. 

We  assume  that  the  material  is  isotropic.  So,  in  the  setting  of  plane  deformation, 
the  stress-strain  relation  is  given  by 

a  =  Atr  (e)I  +  2 ye 
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V  \  dy  “T  dx  J 
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where 


e=I(Vu+(Vu)T) 


dui 

dx 

1  ( d%i\  ,  du2  A 

2  \  dy  '  dx  J 


is  the  linear  strain,  I  is  the  2x2  identity  tensor, 


1 

2 


( 8u\  .  du2  \ 
\  dy  '  dx  J 
duy, 
dy 


E 

l"~  2(1  + u) 


and 


A  = 


Ev 

(1  +  v)(l  —  2v) 


(1.5) 


(1.6) 


are  the  Lame  constants,  E  is  Young’s  modulus,  and  v  is  Poisson’s  ratio. 

In  modeling  a  two-phase  elastic  material,  we  assume  that  all  the  material  parameters 
g,  A,  E,  and  v  are  piecewise  constants.  In  particular,  we  assume  that  the  shear  modulus 
and  Poisson’s  ratio  are  given,  respectively,  by 


f  g+  in  fl+ 
\  g~  in  f 


and 


j  in 
\  v~  in  ’ 


where  all  g+ ,  g  ,  z'+,  v  are  positive  constants.  Usually,  Poisson’s  ratio  v+,v  € 
(0,0.5)  (cf.  [1,  11]). 

We  let  /  =  —Fi/(n  +  A)  and  g  =  —F^Kg  +  A),  and  equations  (1.1)  and  (1.3)  can 
be  written  as 


xd2ui  „  .  d2ui  d2u2  . 

2(!  -  v) lu  +  U  “  2u)+E7T  +  u-l:  = 


(1  -  2u) 


dx2 

d2U2 


dy 2 
d2U2 


dx 2  +  2(1  ^  dy 2 


+ 


dxdy 

d2ui 

dxdy 


2g 


1  -  2v 


=  g(x,y), 


du^. 


,  dui  du2\  ( du\ 

1  “  +  U~E~  ni  +  f*  I  n2 

dx  dy  J  \  dy  dx 


dui  du2\  2g 


dui  du2  , 

■YW (1  ~  n2 


=  i>i 


(1.7) 

(1.8) 

(1.9) 

(1.10) 


where  [•]  is  defined  as  the  jump  of  the  quantity  between  the  outside  and  the  inside  of 
the  interface. 

There  exist  several  numerical  methods  for  solving  general  elasticity  problems  that 
do  not  involve  interfaces.  Among  them,  the  finite  element  method  and  the  boundary 
integral  or  boundary  element  method  appear  to  be  very  successful,  cf.  e.g.,  [3,  19,  23] 
and  the  references  therein.  However,  in  treating  moving  interface  problems,  the  use 
of  fixed  Cartesian  grids  often  shows  advantages  in  practical  computations  [12].  It  is 
therefore  desirable  to  develop  new,  efficient  methods  based  on  finite  difference  schemes 
on  fixed  Cartesian  grids  for  our  underlying  elasticity  problems  with  interfaces.  This  is 
our  primary  goal  of  the  work.  Our  idea  is  to  use  the  immersed  interface  method  [22, 
15,  16]  to  derive  a  finite  difference  scheme  for  the  elasticity  problem  with  an  interface. 
This  is  a  natural  strategy,  since  the  geometrical  complexity  of  the  problem  is  local. 

The  curved  interface  in  the  underlying  problem  brings  up  several  substantial  diffi¬ 
culties  in  the  development  and  analysis  of  numerical  schemes. 
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•  Discretization.  With  a  uniform  Cartesian  grid,  the  interface  is  typically  not 
aligned  with  the  grid  but  rather  cuts  between  grid  points.  Thus,  for  grid  points 
near  the  interface,  the  stencil  of  a  standard  finite  difference  scheme  will  contain 
points  from  both  sides  of  the  interface.  Due  to  the  non-smoothness  of  u\  and 
U‘2-  differentiating  u\  and  u2  across  the  interface  using  standard  finite  difference 
schemes  will  not  produce  accurate  approximations  to  the  derivatives  of  u\  and  u2. 
The  cross  derivative  terms  in  the  differential  equations  need  special  treatments 
in  the  discretization  on  the  grid  points  near  or  on  the  interface. 

•  Solving  the  system  of  discrete  linear  equations.  Because  of  the  presence  of 
interfaces  and  non-smoothness  in  the  solution,  the  system  of  discrete  equations 
is  no  longer  symmetric  or  positive  definite.  The  structure  of  the  system  makes  it 
hard  to  find  an  efficient  solver. 

•  Error  analysis.  Due  to  the  complexity  of  interfaces  and  the  non-smoothness  in 
the  solutions,  it  is  difficult  to  perform  convergence  analysis  in  the  conventional 
way. 

Our  main  contribution  of  this  paper  is  the  development  of  the  immersed  interface 
method  (IIM)  to  the  underlying  elasticity  problems  with  interfaces.  The  key  in  our 
method  is  to  utilize  the  local  coordinate  transformation  to  carefully  analyze  the  rela¬ 
tions  of  quantities  from  one  side  to  the  another.  Such  relations  shall  lead  us  to  find 
accurate  finite  difference  schemes  for  grid  points  near  or  on  the  interface. 

The  paper  is  organized  as  follows.  In  Section  2,  we  describe  our  algorithm.  In 
Section  3,  we  introduce  local  coordinate  transformation  that  is  essential  to  the  devel¬ 
opment  of  our  method.  In  Section  4,  we  describe  interface  relations.  In  Section  5,  we 
derive  the  finite  difference  scheme  for  irregular  grid  points.  The  linear  solvers  are  also 
discussed.  In  Section  6,  we  present  our  numerical  results.  Finally,  in  Section  7,  we 
draw  conclusions. 

2  Description  of  the  Algorithm 

For  simplicity,  we  assume  that  the  domain  D  is  a  square:  D  =  (a,  b)  x  (c,  d )  with 
d  —  c  —  b  —  a.  Let  n  >  1  be  an  integer  and  h  =  (b  —  a)  jn  =  (d  —  c)  jn.  Let 

Xi  =  a  +  ih,  y3  =  c  +  jh ,  i,  j  =  0,  ■  ■  ■  ,  n. 

We  wish  to  solve  the  problem  using  a  finite  difference  method  on  the  uniform  Cartesian 
grid.  Our  result  will  be  a  finite  difference  scheme  of  the  form 

E  ak{Ui)i+ik,j+jk  +  'E  Pk(u2)i+ik,j+jk  = 

k  k 

Elk{U\)i+ikj+jk  +  E  ^k{U'2)i+ik,j+jk  = 
k  k 

at  any  grid  point  (xi,yj),  where  and  {U2)ij  approximate  ui(xl.yJ)  and  u2(xi,yj), 

respectively,  ftJ  =  f{xu  y:j),  g%3  =  (j{xl.y:j).  and  all  ak.  /3k,  7*,  and  5k  are  undetermined 
coefficients.  Each  sum  over  k  in  (2.1)  involves  only  finite  number  of  grid  points  that 

are  centered  at  (xi,yj)  (at  most  nine  grid  points  are  involved  in  our  algorithm),  all 

ik.jk  €  {—1, 0, 1}.  The  coefficients  ak.  j3k,  7^,  and  5k  and  the  indices  ik,  jk  will  depend 


fa  +  ct, 


9ij  +  Cfj, 


(2.1) 
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on  We  omit  the  dependency  for  simplicity.  By  finding  these  coefficients  properly, 

we  can  obtain  a  second-order  accurate  finite  difference  scheme. 


2.1  Classification  of  grid  points 

Before  we  explain  the  finite  difference  scheme,  we  classify  all  grid  points  into  two 
categories:  regular  and  irregular.  We  say  a  grid  point  ( Xi,yj )  is  regular  (Figure  2.1  (a)) 
if  the  interface  T  does  not  cut  through  any  points  in  the  standard  nine  point  stencil 
centered  at  (xj  ,yj).  We  say  a  grid  point  is  irregular,  if  it  is  not  regular.  An  irregular 
grid  point  is  further  classified  as  type  I,  if  the  interface  crosses  the  five  point  stencil 
centered  at  this  point  (Figure  2.1  (b)),  and  type  II,  if  otherwise  (Figure  2.1  (c)  and 

(d))- 


Figure  2.1:  Classification  of  grid  points. 


At  a  regular  grid  point  (Figure  2.1  (a)),  we  use  the  standard  central  finite  difference 
scheme.  At  an  irregular  grid  point  (see  Figure  2.1  (b)-(d)),  we  derive  a  finite  difference 
scheme  according  to  how  the  interface  cuts  through  the  five-point  stencil. 
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2.2  Discretization  at  regular  grid  points 

At  a  regular  grid  point  ( Xi,yj ),  we  have  the  following  approximations  of  the  second- 
order  partial  derivatives  for  a  given  smooth  function  u 

=  "fe.w-i)-Mxi,w)  +  a(xi,w+i)+0(/i2)| 

^  _  u{xi+1,yj+1)  +  u{xi-i,yj-i)  -  u(xi-i,yj+i)  -  u{xi+i,yj-i)  , 

uxy\xi->yj)  ~  4:h?  '  \J\n  J, 

cf.  Figure  2.2.  These  lead  to  the  following  discretization  of  (1.7)  and  (1.8): 

(  (U^j  +  (C7i)i+i j)  +  l~^  (  (Uikj-i  +  (Cfi)iJ+i) 

- ^2 - +  4^2  (  (^i+lj+l  + 

—  (U2)  i—  lj'+l  —  ( ^2)i+l, j—1  )  =  fiji 

(2.2) 

—^2 —  (  (^2)i-l, j  +  (^2)i+l,j)  d - ^2 -  (  (^2)ij-l  +  (^2)i,j+l) 

- J^2 - +  4^2  (  (^l)i+lj+l  +  (^l)i-lj-l 

—  (C7i)i-ij+i  —  (Li)j-i-ij— 1  ^  =  gij , 

where  {U)ij  stands  for  the  approximation  of  u(xi,yj ),  and  in  this  case  C\3  =  Cfj  =  0, 
cf.  (2.1). 


Figure  2.2:  Discretization  at  a  regular  grid  point. 


2.3  Discretization  at  type  I  irregular  grid  points 

Following  [15]  and  [16],  we  use  a  nine-point  stencil  for  u\  and  it 2  in  both  equations  of 
(2.1)  at  such  a  grid  point  (xi,yj)  (cf.  Figure  2.1  (b)).  To  determine  the  coefficients 
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of  the  finite  difference  equations,  we  first  choose  a  point  (x*,y* )  on  the  interface  T 
that  is  near  ( Xi,yj ).  We  replace  (Ui)i+ikj+jh  and  (U2)i+ikj+jh  with  the  exact  solution 
ui{xi+ik,yj+jk)  and  U2{xi+ik,yj+jk)  in  (2.1)  and  use  the  Taylor  expansion  at  (x* .  y* ) 
for  each  term  in  order  to  set  up  a  system  of  equations  for  the  coefficients.  Since  the 
solution  from  both  sides  is  involved,  we  use  the  superscript  —  and  +  to  denote  the 
limiting  values  of  a  function  from  one  side  or  the  other.  For  example,  if  the  point 
(xi,yj)  is  located  inside  the  interface,  then  we  can  expand  a{xl,y.j)  to  get 

u{xi,yj)  =  u~  +  u~(xi  -  x*)  +  u~  (yj  -  y*)  +  ^uxx{xi  -  x*)2 

+  \uyy(yj  -  y*j )2  +  u~y{xi  -  x*)(yj  -y*j)  +  0{h 3). 

If  we  do  such  expansion  at  each  grid  point  used  in  the  finite  difference  equations  in 
(2.1),  then  the  local  truncation  errors  Tf  and  T-j  (for  the  first  and  second  equation 
respectively)  can  be  expressed  as  a  linear  combination  of  all  the  values  uf,  ufx,  ufy, 

ulxxi  '  '  '  5  u2 yyi  an(^  u2xy 

Next,  we  express  all  the  values  from  +  side,  uf ,  ufx.  ufy.  •  •  • ,  ufyy,  ufxy.  in  terms 
of  the  values  on  the  —  side,  wj”,  px,  u~[y,  ■  ■  ■ ,  U2XX,  p, yy,  up.  To  do  so’  we  need  to 
use  the  interface  conditions, 

uf  =  p,  uf=u  2,  (2.3) 

1  (  (!  -  v+)(ui)t  +  v+(u-2)f )  ni  +  p  (  (ki)+  +  (ti2)+)  n2 

=  (  (!  -  *0(«i)*  +  ”~(u2)y)  ni 

+  P  (  («i)y  +  («2)x)  n2  +  </>,  (2-4) 

/^+  (  («i)»  +  («2)j)  «i  +  j ~rp+  (  v+(ui)t  +  (1  -  ^+)(w2)y  )  n2 

=  (  («i)»  +  (“2)*)  m  +  1  2_^2u_  ( 

+  (1  -  W)(rt2)“)  n2  +  -0.  (2.5) 

To  obtain  all  the  coefficients  a*,,  /3k,  7 and  <5*  in  the  finite  difference  equations  (2.1), 
more  interface  relations  are  needed.  Differentiating  the  above  equations  and  manipu¬ 
lating  the  results  allow  us  to  perform  the  desired  elimination,  as  detailed  in  Section  4. 
In  order  to  do  so,  it  is  very  convenient  to  first  perform  a  local  coordinate  transformation 
in  the  normal  direction  £,  and  the  tangential  direction  r]  of  T  at  (x*.  y*  ). 

Once  the  local  truncation  errors  Tf  and  Tf  are  expressed  as  a  linear  combination 
of  the  values  from  just  one  side,  say  tij",  ulx,  ujj  ,  pxx,  ujj  ,  pxy,  u%  ,  up  u%y,  up, 
u2yy,  u2xy,  we  must  require  that  the  coefficient  of  each  of  these  terms  to  vanish  in  order 
to  match  the  partial  differential  equation  up  to  second  order  derivative  terms.  This 
gives  us  a  system  of  twelve  linear  equations  from  the  first  and  second  finite  difference 
equations  respectively.  Nine-point  stencil  will  be  used  to  obtain  a  solvable  system,  see 
Section  5  for  all  these  derivations. 
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2.4  Discretization  at  type  II  irregular  grid  points 

For  a  type  II  irregular  grid  point  (cf.  Figure  2.1  (c)  and  (d)),  we  use 

•  five-point  stencil  for  u\  and  seven-point  or  four-point  stencil  for  the  cross  deriva¬ 
tive  of  u-2  in  the  first  equation,  and 


•  five-point  stencil  for  u2  and  seven-point  or  four-point  stencil  for  the  cross  deriva¬ 
tive  of  u\  in  the  second  equation. 

Therefore  equations  (1.7)  and  (1.8)  have  the  following  form: 


2(1  —  u) 


(  {Ui )i-u  +  (C7i)i+ij)  +  (  (Cfi)ij-i  +  (Ui)i,j+i) 

2{*~2Au\ui)i,j  + 


1  -  2v 


K2 


(  (U2)^ij  +  (U2)i+ij) 

-  h-2 — "(^ )i,j + y^7fc(^i)i+u 


d+3k  —  fiji 

2(1  —  v) 


(  (U2 )iJ_1  +  (U2)iJ+ 1) 


d+jk  ~  9ij- 


(2.6) 


(2.7) 


To  determine  the  coefficients  f3k  and  7 k  in  (2.6)  and  (2.7),  we  use  four-point  stencil. 
For  the  grid  points  shown  in  Figures  2.1  (a)  and  (b),  we  can  apply  the  following 
expansion  for  a  smooth  function  u  to  derive  our  finite  difference  scheme 


^ xy 


ui1)  —  ui2)  —  U +  U 

h2 

U ^  —  U —  U ^  +  U ^ 

K2 

U^  —  U —  U ^  +  U 

h2 

UW  -  u («)  - 
h 2 


+  0{h) 

(2.8) 

+  0{h ) 

(2.9) 

+  0(h ) 

(2.10) 

+  o(h), 

(2.11) 

where  i/W  denotes  u  value  at  the  point  numbered  i  in  Figure  2.1.  Which  formula  from 
(2.8)  (2.11)  to  use  depends  on  how  the  interface  crosses  the  nine  -point  stencil.  In 

the  specific  case  in  Figure  2.1  (d),  we  use  the  expansion  (2.10)  or  (2.11),  where  points 
numbered  (4)  to  (9)  are  in  the  same  side  of  the  interface.  However,  when  one  grid  point 
is  on  the  other  side  of  the  other  eight  grid  points  in  the  nine-point  stencil  as  shown  in 
Figures  2.3  (a)  and  (b),  we  apply  the  following  formula  for  a  smooth  function  u.  In 
this  case  we  have 


Uxy 


Ui2)  —  u +  u —  2  u^  +  u —  u +  u 
~2h2 

^i1)  _  ui2)  —  ui4)  +  2 ui5)  —  u^)  —  ui8)  +  u )9) 

2h? 


+  0{h2) 

(2.12) 

+  0{h2). 

(2.13) 

Whether  we  use  (2.12)  or  (2.13)  depends  on  the  geometry.  We  use,  for  example,  the 
formula  (2.12)  in  Figure  2.3  (a),  and  (2.13)  in  Figure  2.3  (b). 
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(a) 


(b) 


Figure  2.3:  An  illustration  of  the  geometry  of  type  II  irregular  grid  points. 


3  Transformations  of  Local  Coordinates 


For  each  type  I  irregular  point  (xi,yj),  we  need  to  find  a  point  (x*.y*)  on  the  interface. 
We  usually  take  this  point  as  the  projection  of  (xi,yj)  on  the  interface  if  the  interface 
is  smooth  at  this  point.  Otherwise  we  can  take  any  point  on  the  interface  where  the 
interface  is  smooth. 

After  (x*.y*)  is  selected,  we  apply  the  local  coordinate  transformation  at  this  point. 
Let  #  be  the  angle  between  the  x— axis  and  the  normal  direction,  pointing  in  the 
direction  of  the  +  side,  cf.  Figure  3.1.  The  transformation  is  defined  as  follows: 

\  t  =  (x  ~  xt)  cos($)  +  {y  —  yj)  sin(0), 

<  (3.1 

(  r?  =  -(x  -  x*)  sin(#)  +  (y  -  y*)  cos(#). 

Under  this  local  coordinate  transformation,  the  governing  equations  in  (1.7)  and  (1.8) 
become 

(cos2  #  +  1  —  2i')'iii££  —  2  sin  #  cos  9u\^v  +  (sin2  9  +  1  —  2v)uim 
+  sin  #  cos  0U2££  +  (cos2#  —  sin2#)^^  —  sin#cos#U2W  =  /,  (3.2) 

(sin2  9  +  1  —  2v)v,2££  +  2  sin  9  cos  9u2^v  +  (cos2  9  +  1  —  2u)u2 nn 
+  sin#cos#ui££  +  (cos2#  —  sin2#)tti^  —  sin  #  cos  9u\vv  =  g1  (3.3) 

where  /  =  f(x,y)  =  /(£,??),  g  =  g{x,y)  =  g{^y)  for  simplicity. 


Figure  3.1:  The  local  coordinate  transformation. 
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4  Interface  Relations 

We  consider  a  fixed  point  (x* ,  y* )  and  define  a  new  £-r]  coordinate  system  based  on 
the  directions  normal  and  tangential  to  T  at  this  point  using  the  formulas  (3.1).  In 
a  neighborhood  of  this  point,  the  interface  can  be  parameterized  as  £  =  %(ry),  y  =  y. 
Notice  that  x(0)  =  0  and  x'(0)  =  0  provided  that  the  interface  is  smooth  at  {x*,y*). 
We  need  to  express  all  the  quantities  with  the  superscript  +  in  terms  of  those  quantities 
with  the  superscript  — . 

First,  using  the  same  approach  as  that  in  [15]  and  [16],  we  easily  get  the  following 


jump  conditions  for  u\  and  t<2, 

ut(x(v),v)  lu=o  =  «r(x(»y)>»y)lu= o,  (4.i) 

«2 (x(v),v)\n=o  =  ^2  (x(v),v)\v=o*  (4-2) 

[«if]x"  +  =  0,  Mx"  +  [u2nn]  =  0.  (4.3) 

The  jump  (4.3)  can  be  rewritten  as 

uivv  =  uirin  +  -  x”'>4v  (4-4) 

U2VV  =  u2 nn  +  x"u2i  ~  x"«2f  (4-5) 


Notice  that,  in  the  local  coordinate  system,  the  unit  normal  at  any  point  (x(r;),  y) 
on  the  interface  is 


Therefore,  at  (x|,yj),  or  (0,0)  in  the  local  coordinate  system,  we  have 
cr-  i  _  X'(V)X”(V) 

nlr]\ri=0  ~  ssos  5  ’ 

(1  +  {x'iv))2)*  n= o 

T/=0 

Recalling  that  x(0)  =  0  and  x'(0)  =  0,  we  have,  at  (x^,y*),  or  (0,0)  in  the  local 
coordinate  system,  that 

m  =  1,  «2  =  0,  (4.7) 

nTij  =  0,  =  — X"(0).  (4.8) 

In  the  x-y  coordinate  system,  we  have  n  =  (ni,ri2)T.  The  relation  between  (ni,ri2)T 
and  (ffi,  njj7  is  determined  by 

m  =  ni*  cos  6  —  U2  *  sin  0, 
n2  =  n{*  sin  0  +  *  cos  6. 

Hence,  by  (4.7)  and  (4.8),  we  obtain 


-x''{v)  +  ( x'{v))2x”{v ) 

4  +  Ix'iv))2  v=0  (1  +  (x'(v))2) 1 


n\  =  cos  6 , 


n-2  =  sin  9, 


(4.9) 
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niv\ri=o  =  nirj  1^=0  cos  #  -  n2r]  1^=0  sin  9  =  x"(0)  sin  #,  (4.10) 

n2rt\y=o  =  0^=0  sin#  +  02^=0  cos#  =  — xw(0)  cos  #.  (4.11) 

These  relations  will  be  used  to  derive  the  interface  relations  for  u and  . 

The  following  result  is  required  to  find  the  interface  relations  so  every  quantity  from 
ui^  ’  u%rr  u2£ri ■  ui&  an<4  can  13(3  expressed  in  terms  of  ,  u^x,  u~[y,  u~[xx,  u^yy, 
uixv’  u2j  u2x->  u2v>  u2xx-  U2'wu ■  an<4  u2xv  ^s  proof  is  trivial,  and  is  omitted. 


Lemma  4.1.  If  Det  = 


o  n  ai2 
021  0,22 


0,  then  the  system  of  linear  equations 


j  aux  +  ai2y  =  buz i  +  612^2  H - h  blnzn 

\  a2\x  +  a22U  =  b2\Z\  +  b22z2  H - h  b2nzn 

has  the  unique  solution 


1  ( 

bn 

«12 

Z\  + 

#12 

012 

%2  +  • 

•  + 

bln 

012 

x=DFty 

#21 

022 

#22 

022 

#2n 

022 

1  ( 

an 

#11 

Zl  + 

an 

#12 

Z2  +  • 

•  + 

Oil 

bln 

y  =  Tmy 

<*21 

#21 

<*21 

O* 

to 

to 

021 

#2n 

4.1  Interface  relations  for  and 

Rewrite  the  interface  conditions  (2.4)  and  (2.5)  as  follow: 

ni(a+u+x  +  /3+a+)  +  n2y+{u^y  +  11+) 

=  m  (a~uix  +  p-u^y)  +  n2n~  (u^y  +  «JX)  +  <p{x,y), 

mn+{uiy  +  u^x)  +  n2(f3+Uix  +  a+u%) 

=  ni/i“(«ty  +  ujx)  +  n2(f3~u^x  +  a~  u^y)  +  ip{x,  y ), 

where  n  =  ( n\,n2)T  =  {ni{£,r]),n2(£,ri))T  is  the  unit  normal  to  the  interface  and 


Q+  =  r^F(1  - ^  ^  = 


2y+ 

—  2u+ ' 
2  n~ 


Notice  that 


-7—  =  u^x  +  Uivr]x  =  cos(#)tq {  -  sin(0)ulv, 
Qx 

—  =  ui^y  +  uivr)y  =  sin(#)tti?  +  cos  (#)«!,,, 

^  =  u2^x  +  u2r]rix  =  cos(#)a2£  -  sin {6)u2rn 
Qx 

=  u2  +  u2vr]y  =  sin(#)u2£  +  cos  {0)u2v. 


Thus,  in  the  local  £-rj  coordinate  system,  we  have 

(nia+  cos  9  +  n2y+  sin#)ri^  +  (— nia+  sin#  +  n2y+  cos#)rt^ 


12 


Xingzhou  Yang,  Bo  Li,  and  Zhilin  Li 


+  ( n\f3+  sin  0  +  ri2/r+  cos  9)u^  +  {n\/3+  cos  0  —  ri2/r+  sin  0)it^ 

=  ( n\a~  cos 0  +  ri2/r_  sin 9)u^  +  (— n\a~  sin0  +  n 2/r_  cos 0)u^rj 
+  ( ni/3~  sin  0  +  n 2H~  cos  9)u ^  +  (rii/3-  cos  9  —  n 2/r_  sin  9)u ^ 

+<!>(£, v),  (4-12) 

(/i+rii  sin0  +  ri2/3+  cos0)n^  +  (ni/r+  cos0  +  n2«+  sin0)n^ 

+  (ni//+  cos  0  —  ri2 /3+  sin  0)0^  +  (— /r+ni  sin0  +  n2<a;+  cos  9)u2V 
=  (n~ni  sin0  +  ri2/3_  cos  9)u^  +  (rii/i-  cos  9  +  n2a~  sin  9)u^ 

+  (nin~  cos  9  —  U2/3~  sin 9)u^v  +  sin0  +  r^a-  cos  0)0^ 

+mv)-  (4.13) 

Using  the  fact  that  ni  =  cos  0,  ri2  =  sin  9  at  (0, 0)  in  the  local  coordinate  system 

and  the  jump  conditions  (4.1)  and  (4.2),  we  can  solve  for  and  u.^  from  equations 

(4.12)  and  (4.13)  in  terms  of  ojj,  •  •  •  ,u^  .  -*-11  feet,  using  (4.9)  we  can  rewrite  the 
equations  (4.12)  and  (4.13)  as 

(a+  cos2  9  +  n+  sin2  9)u +  (/ 3+  +  /r+)  sin^cos^n^ 

=  ( a~  cos2  9  +  fi~  sin2  0)u  +  (—a-  +  n~)  sin  9  cos  00^ 

—  (— a+  +  /i+)  sin  9  cos  00^  —  (f3+  cos 2  9  —  n+  sin 2  9) 

+  ( f3~  +  n~ )  sin  9  cos  9u ^  +  (/ 3~  cos2  9  —  /j~  sin2  9)^^ 

+4>{^V),  (4-14) 

(n+  +  (3+ )  sin  9  cos  9u^  +  (/r+  cos2  9  +  a+  sin2  9)u ^ 

=  ( fi~  +  j3~)  sin 9  cos  9u^  +  (/j~  cos2  9  +  a~  sin2  9)u ^ 

—  (n+  cos2  9  —  /3+  sin2  0)oj^  —  (— /r+  +  a+)  sin  9  cos  00^ 

+  (/r_  cos2  9  —  (3~  sin2  9)u +  (— /r_  +  a-)  sin  0  cos  0 

(4-15) 

Since  it]1"  =  u ^  and  =  u^,  the  right  hand  sides  of  (4.14)  and  (4.15)  become 

b\  =  (a-  cos2  9  +  n~  sin2  9)u^ 

+  [(a+  —  a~)  —  (fi+  —  n~)\  sin0  cos  9u^v  +  (/ 3~  +  n~)  sin  0  cos  00^ 

+  [{H+  —  n~)  sin2  9  —  (/3+  -  f3~)  cos2  0]^  +  <!>(£,  r)), 

02  =  (/r_  + /3_)  sin0cos0uj~£ 

+  [(/3+  —  (3~)  sin2  9  —  (/r+  —  n~)  cos2  9]u{v  +  (ji~  cos2  0  +  a~  sin2  9)u^ 

+  [(/r+  -  n~)  -  (a+  -  a-)]  sin  0  cos  07/^  +  ^(£,0)- 

(4.16) 

Let 

au  =  cos20a+  +  sin20/r+,  ai2  =  {/3+  +  n+)  sin  0  cos  0, 

021  =  (/r+  +  /3+ )  sin  0  cos  0,  022  =  0+  cos2  0  +  a+  sin2  0. 
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The  determinant  of  the  coefficient  matrix  for  the  linear  system  (4.14)  and  (4.15)  is 


an 

012 

021 

022 

cos2  9a+  +  sin2  #/r+ 
+  /3+)  sin  9  cos  9 


(/3+  +  n+)  sin  9  cos  9 
H+  cos2  9  +  a+  sin2  9 


=  or  iv 


Since 


a+ n+  = 


2(#+)2 

1  -  2u+ 


(1 


-  ^+)  f  o, 


(4.17) 


the  system  (4.14)  and  (4.15)  has  the  unique  solution 


,+  _ 

1 

bi 

012 

'If  ~ 

a+/i+ 

b2 

022 

,+  _ 

1 

an 

bl 

'2?  “ 

a+/i+ 

021 

b2 

1 

a+n+ 

1 

a+/i+ 


(022&1  —  Ol2&2), 


(an&2  —  <*2i&i  • 


(4.18) 

(4.19) 


Consequently,  the  values  u ^  and  can  be  expressed  in  terms  of  .  u^,  ■  ■  ■ ,  . 

Moreover,  plugging  (4.18)  and  (4.19)  into  (4.4)  and  (4.5),  we  can  also  express  and 
u2m  in  terms  of  those  values  in  the  —  side. 


4.2  Interface  relations  for  and 

To  find  the  interface  relations  for  and  u2^r)-,  we  differentiate  (4.12)  and  (4.13)  with 
respect  to  r).  After  some  calculations,  we  obtain  the  following  system  of  linear  equations 
about  and  u2^ 


M\ui +  Ai2'u+?r(  =  B  i, 

+  ^22^2^7)  =  -®2, 


(4.20) 


where  the  coefficients  An,  A12,  A21,  A22  are  given  by 


An  =  niQ.+  cos  9  +  n2fi+  sin#,  A12  =  ni/3+  sin#  +  712/1 +  cos  9, 

A-21  =  nin+  sin  9  +  ri2/3+  cos  #,  A22  =  n\ /i+  cos  9  +  ri2a+  sin  9. 

The  terms  B\  and  B-2  involve  u^.  u2^,  u2r)V,  and  u7  ,  u^,  ■■■,  .  However, 

since  we  have  already  obtained  that  u\^.  u2^ ,  u\rrq,  and  u2rm  can  be  expressed  as  a 
linear  combination  of  u~[ .  u^,  ■  •  • ,  ,  the  values  B\  and  B2  can  also  be  expressed  in 

terms  of  these  quantities  as  well.  So  we  set 

Hi  =  b\u7  +  b\u^  +  b\uUj  +  b\u^  +  bluUjrj  +  b\u^ 

+b\u2  +  b\u~^  +  b\u2r)  +  b\0u~ K  +  blnu2m  +  b\2u2^  +  6j, 

B2  =  bfu^  +  +  bjuUl  +  +  bju^ 

+b27u2  +  +  blu2r}  +  +  hliu2m  +  &?2^  +  bl 

Since  71 1  =  cos  9  and  712  =  sin#,  the  determinant  of  the  coefficient  matrix  of  the  system 
(4.20)  is 


An 

A12 

cos2  9a+  +  sin2  #/i+ 

(/ 3+  +  /i+)  sin  #  cos  # 

A21 

A22 

(n+  +  (3+)  sin# cos# 

/i+  cos2  #  +  a+  sin2  # 

DET  = 


a+n+  /  0. 
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Thus,  from  Lemma  4.1,  the  solution 
1 


u 


+  _ 


U 


DET 

1 

+ 


b\ 

b\  A 


+ 


+ 


+ 


DET 

1 

DET 

1 

DET 

1 

DET 
1 


A\2 
l22 
b\  A12 
b\  A22 

A12 
A22 

A12 
bjo  A22 
b\2  A\2 
‘-'22 


Ux  + 


to  (4.20)  is 
1 


DET 


b\ 


q  A 


A12 

22 


UlZ  +  DET 


% 

% 


A12 

A22 


u 


'1 V 


b} 

b2r 

b\ 0 


u^  + 

+ 


1 


u. 


DET 

1 

DET 
1 


b\ 

^4l2 

1 

7  /  I 

q 

^4i2 

q 

A22 

W'  DET 

q 

A22 

bl 

b2s 


DET 
1 


+ 


+ 


+ 


+ 


DET 

1 

DET 

1 

DET 
1 

DET 


bl  2  A 

An  b\ 
A21 
An 
A21 


+  DET 
1 

+  DET 


A 12 
A22 

A  l2 
bh  A22 
bl  A12 

q  A22 


U'2t  +  DET 


bl  A, 


bl 


'9  -^12 

X22 


u 


'2  V 


b\  1 


u 


2m 


(4.21) 


bl 


Ui  + 


DET 


An  bl 
A21  b2 


+  DET 


An  bl 

A21  q 


’hr) 


bl 

bl 


+  DET 


An 

q 

1 

An 

q 

A21 

q 

lr>n  ^  DET 

A21 

q 

An  b\ 


A21  q 

All  frio 
A2i  b\  0 

A11  b\ 2 


Uo  + 


1 


DET 
1 


An 

q 

1 

1 1  1 

An 

q 

A21 

q 

2«  +  DET 

A21 

q 

u 


2v 


+  DET 


^21 


bn 


"2 +  DET 


1 


An  b\  1 
^21  bh 
An  b\ 

^21  q 


u. 


'2  m 


(4.22) 


These  are  the  desired  relations. 


4.3  Interface  relations  for  and 

From  the  differential  equations  (3.2)  and  (3.3),  we  have 

(cos26»  +  l-2i/+)tt+?-2  sin  9  cos  +  (sin20  +  1  -2u+)u\ 

+  sin  0  cos  9u ^  +  (cos2  0  —  sin2  9) —  sin  9  cos  60  ^  —  f+ 

=  (cos 2  9  +  1  —  2v~)u ^  —  2  sin#  cos  60^  +  (sin2#  +  1  —  2u~)uirjrj 

+  sin  9  cos  9u^  +  (cos2  9  —  sin2  0)uj^  —  sin  9  cos  9^^  —  f~,  (4.23) 

(sin2  9  +  1  —  2v+)v,2££  +  2  sin  9  cos  9u^v  +  (cos2  9  +  1  —  2u+)u2m 
+  sin  9  cos  9u^  +  (cos2  9  —  sin2  9)u^r)  —  sin  9  cos  9u^m  —  g+ 

=  (sin2#  +  1  —  2u~)u~ii  +  2  sin 9  cos  Ou^  +  (cos 2  9  +  1  —  2i2~)u2T]r} 

+  sin  9  cos  9u^  +  (cos2  9  —  sin2  9)u^  —  sin  0  cos  9u^m  —  g~ .  (4.24) 

We  need  to  solve  the  above  equations  for  and  Observe  that  the  determi¬ 
nant  of  the  coefficient  matrix  of  the  above  system  is 


1  —  2u+  +  cos2  9  sin  9  cos  9 

sin  9  cos  9  1  —  2v+  +  sin2  9 


2(1-2i/+)(1-i/+)  ^0, 
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By  (4.4),  (4.5),  (4.21),  and  (4.22),  we  can  solve  (4.23)  and  (4.24)  for  and 
which  are  expressed  in  terms  of  tt^,  u^,  •  •  • , 

5  Derivation  of  the  Finite  Difference  Scheme  at  Type  I 
Irregular  Grid  Points 

We  are  now  ready  to  derive  the  finite  difference  scheme  for  irregular  grid  points  of  type 
I.  We  first  rewrite  the  finite  difference  scheme  (2.1)  as  follows: 

{9  9 

X]  ak{Ul)i+ikj+jk  +  X]  f3k{U'2)i+ik ,j+jk  =  f{xi,yj)  +  C}p 

kt  *9=1  (5-1) 

X/  lk{Ul)i+ik ,j+jk  +  X]  $k{U2 )i+ik,j+jk  =  g{Xi,Vj)  +  Cfj. 

k=l  k=l 

We  use  the  undetermined  coefficient  method  to  determine  all  the  coefficients  ak ,  /3/ c, 
7fc,  and  Sk. 

Denote  the  £-7/  coordinates  of  the  nine  grid  points  in  the  finite  difference  stencil 
(see  Figure  5.1.) 

{xi-i,yj-i),  (zj-1,7/7),  {xi-i,yj+1),  (xi,yj- 1),  (zj,7/j), 

{xi, yj+ 1),  (2^+1, yj— 1),  (2^+1,  yj):  (2^+ 1, y^+i) 
as 

(£l,T7l),  (6:^2),  (6,773),  (£4,774),  (£5,775),  (£6,  776),  (£7,777),  (£s,  77s),  (£9,779)- 


Figure  5.1:  The  labels  of  the  nine  grid  points. 

The  local  truncation  errors  at  a  grid  point  ( Xj^yj )  are  defined  as 
9  9 

=  X]  afcKl  ,  %+jJ  +  X]  ^*«2(a:i+i* ,  ?7j+jJ  -  f{xi,yj)  ~  C}j 

k= 1  /c=l 

9  9 

=  ^ayfc77l(£/c,7?/c)  +  ^PkU2  {(,k,Vk)  ~  fithVj)  ~  Qj, 

k= 1  /c=l 


(5.2) 
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Tfj  =  ^lkui{xi+ik,yj+jk)  +  ^sku2{xi+ik,yj+jk) -g(xi,yj)  -  Cg 

k  i  k  1 

9  9 

=  +  ^2skU2{^k,Vk) -gitiiVj)  -  Cfj.  (5.3) 


fc=a 


jfc=i 


We  expand  u\ {€k,Vk)  and  u^i^k-  Vk)  about  (0, 0)  in  Taylor  series  in  the  local  coordinate 
system  from  each  side  of  the  interface  to  get 


«i(&»  Vk)  =  uf  +  £ku%  +  Vkutv  +  +  CkJlkU%n  +  \vkutvv 

+0(h3), 

u2^k,  Vk)  =  u2  +  +  ^kVkU^r)  +  rlku2riri 

+0(h3), 


(5.4) 


(5.5) 


where  the  +  or  —  is  chosen  depending  on  whether  (£*,,  r/k)  lies  on  the  +  or  —  side  of  T. 
Carrying  out  such  expansion  for  each  point  involved  in  (5.2)  and  (5.3),  we  obtain  the 
following  expressions  of  and  for  the  local  truncation  errors  as  linear  combinations 
of  the  values  u ^  ^1^’  ^1  nrp  ^2  '  ^2^’  ■  ^2££’  ^2^’  and  u^jyq- 


T1 

u 


rn2 

*3 


aiux  +  a2^i  +  03^  +  +  05^  +  06*^ 

+(27^-^^^  H"  "I-  ^9^1  rjr)  ^10^17777  “I”  "t" 

+  &1«2  +  62^2  +  ^3%^  +  &4«2£  +  &5«2tj  +  ^U2rj  (5-6) 

+b7u^  +  bsu^  +  b9u^vv  +  b10  ufm  +  ftnit^  +  6120^ 

-/,7  -  <71.  +  c?(*), 

CiU^  +  C2tt+  +  c3u^  +  c4n^  +  csu^  +  C6ufv 

+C7«r«  +  C8^i?  +  C9%w  +  ClOUi^  +  +  C12« 

+C?it<2  +  ^2^2  +  ^34*2£  +  ^4^  +  ^5^277  +  ^6^  (5-7) 

+rf7«^f  +  <fe«2e£  +  dgU2vv  +  dl0«aj„  +  rfll«2e„  +  rf12«2£„ 

All  the  coefficients  a*,  6jt,  c*,  and  dk  depend  only  on  the  position  of  the  stencil  relative 
to  the  interface.  In  particular,  they  are  independent  of  the  functions  ui,  112,  /,  and  g. 
If  we  define  the  index  sets  K+  and  K~  by 

K±  =  {k  :  (£&,  Tjk)  is  on  the  ±  side  of  T}, 


then  the  coefficients  a*,  bk ,  Ck,  and  dk  are  given  by 


ai  =  E  otk, 

keK~ 

04  =  E  £k®k, 
keK+ 


02  —  I]  CXk, 
keK+ 

05  =  I]  gkoik, 

keK- 

le2. 


03  =  I]  CfeCtA;, 
keK~ 

06  =  X]  VkOlk i 
keK+ 

1  „2  . 


&7  —  X]  2^kGih'>  X]  2  X)  2 

keK-  keK+  keK- 

Oio  =  I]  \Vkak,  °n  =  X)  ZkWkOlk,  O12  =  E  CkVkOtk, 
k£K+  keK-  k£K+ 


(5.8) 
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bi  =  £  Pm,  b2=  £  Pm, 

m£K~  m£K+ 

bi  =  £  £>mPm,  b§  —  £  IJmPm, 

m£K+  meK~ 

b7=  £  \imPm,  b8=  £  \i2mPm 


£  Pm,  b8  £  £,mPm , 

m£K+  m£K~ 

=  £  VmPm ,  b§  =  £  IJmPm, 

meK~  meK+ 

=  £  himPm,  b9=  £  ^Pm, 

m£K+  m£K~ 


bio  —  £  2  rimPm,  b\  1  —  £  fimVmPm,  ^12  —  £  ^mVmPm, 


meK+ 


meK+ 


Cl  =  £  7 k,  C2=  £  7 k,  C3=  £  67*, 

feeic-  teJf+  keK~ 

C4  =  £  6c7fc;  c5  =  £  %7/c5  c6  =  £  %7fe, 

fceJR'+  kex-  keK+ 

C7  =  £  5^fe7fc:  C8=  £  Klnfk,  Cg  =  £ 

fee/r-  /ce7s:+  k£K~ 

ClO  =  £  5^7^  C11  =  £  tkVklk,  C12  =  £  tkVklk, 

k£K+  k£K~  k£K+ 


(5.10) 


dl  —  £  $m,  d-2  —  £  $m,  C?3  —  £  £,m^m, 

meK~  m<EK+  meK~ 

di  =  £  Cm^m,  d§  =  £  IJm^m,  d%  =  £  T]m$m, 

m£K+  m£K~  m£K+ 

di=  £  \(m^m,  d8=  £  d9  =  £  ^rf2n8m, 

m£K~  m£K+  m£K~ 

dlO  £  2  Vm^m,  d\i  £  £,mVm$m,  di2  £  £,mVm$ri 


m£K+ 


m£K+ 


(5.11) 


Let  U ix ,  Uixx,  Uixy ,  U2x,  U2xx,  U2xy  be  12  dimensional  column  vectors  with  every  co¬ 
ordinate  for  each  vector  corresponding  one  of  the  coefficients  of  wj” ,  t/£,  u~[r).  ujrrn  ■ 

ul £y,  U2  ’  U2^  U2 yi  U2 !££’  U2 riy,  U2(j)  *n  expressions  of  U^,  u2(r)i  u2 jg|- 

We  can  then  rewrite  u+£,  u^.  u^,  in  terms  of  Uix,  Uixx,  Uixy,  U2xi  U2xx, 


U2xy.  Let 


Vec-  = 


Ul  ,  Ul£,  Ul T), 


ul££,  Ulyy,Ul^y,  U2  i  U2 £’  U2yTU2^  U2r)T)iU2£ri 


(5.12) 


From  the  interface  conditions  we  derived  in  Section  4,  we  have  the  following  simple 
expression  in  the  dot  product  form 


u+f:  =  Vec-  ■  Uix  +  wi£0, 

=  Vec-  ■  U2x  +  w2£0, 
^££  =  Vec-  ■  Uixx  +  »i«0, 
«2££  =  Vec-  ■  U2xx  +  w2££0 , 
u ^  =  Vec-  ■  Uixy  +  wifro, 


(5.13) 


-  V CC—  •  U2xy  W2^y0  . 


Using  the  jump  conditions  for  and  it^  in  (4.4)  and  (4.5),  we  can  rewrite  the 
truncation  errors  (5.6)  and  (5.7)  as 


Tij  =  (ai  +  02)%  +  (o3  +  x"aio)n1£  +  (a4  -  x"ai0)n+£ 
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+  (05  +  +  [07  -  (cos2  9  +  1-  2u~)]u^  +  a%u^ 

+  [a9  +  aio  -  (sin2  9  +  1-  2v~)\i+[vv 

+  [an  +  2  sin  9  cos  0]u^rj  +  a12u^n 

+  (&i  +  b2)u2  +  (&3  +  x"bio)u2£  +  (&4  -  x"bio)u^ 

+  (&5  +  b%)u2V  +  (67  -  sin  9  cos  9)u^  +  hu^ 

+  (&9  +  bio  +  sin  9  cos  9)u ^  +  [611  -  (cos2  9  -  sin2  9)]u^v  +  b12u^v 
+  [(cos20  +  1  —  2u~)u ^  —  2  sin  0  cos  9u^rj  +  (sin20  +  1  —  2v~)uivr] 

+  sin  9  cos  9u +  (cos2  9  —  sin2  9)u ^  —  sin  9  cos  9 —  /-] 

-C%  +  0(h) ,  (5.14) 

T?j  =  (ci  +  c2)u^  +  (c3  +  x"cio)u^  +  (c4  -  x!'cio)ui£ 

+  (C5  +  cq)uiv  +  (c7  -  sin#cos0)tt^  +  c8u^ 

+  (09  +  cio  +  sin  9  cos  9)u^vv 

+  [cn  -  (cos2  9  -  sin2  9)]u~ ^  +  c12u^v 

+  (di  +  d2)u2  +  (d3  +  x"dio)u~^  +  (di  -  x"dio)u^ 

+  (d5  +  do)u7,v  +  [dr  -  (sin2  9  +  1  -  2 v~)]u^  +  d8u^ 

+  [d9  +  e?i0  -  (cos2  9  +  1-  2v~)]u2m 
+  (dn  -  2  sin  9  cos  9)u~ ^  +  d12u^v 

+  [(sin20  +  1  —  2 u~)u^  +  2  sin  9  cos  9u^v  +  (cos 2  9  +  1  —  ^~)u2r]ri 
+  sin  9  cos  9u^£  +  (cos2  9  —  sin2  9)u^v  —  sin  9  cos  9 —  g~j\ 

-C%  +  0(h).  (5.15) 

Notice  that  from  the  equations  (3.2)  and  (3.3),  we  have 


(cos20  +  l  —  2v  )u1^  —  2  s'm  9  cos  9u1^rj  +  (sin20  +  1  —  2v  )ulr]r} 

+  sin  9  cos  9u^  +  (cos2  9  —  sin2  9)u ^  —  sin  9  cos  9u2VV  =  f~j , 

(sin2  9  +  1  —  2u~)u2^  +  2  sin 9 cos  Ou^  +  (cos2  9  +  1  —  2^_)t«2W 
+  sin  9  cos  9u ^  +  (cos2  9  —  sin2  9)u^rj  —  sin  9  cos  9u^rjrj  =  g~-. 

From  (5.12)  and  (5.13),  we  can  rewrite  the  truncation  errors  T^,  in  (5.14)  and 
(5.15)  as 


Tlj  =  Vec- 


Li  +  (a4  —  x"aw)Uix  +  a8Uixx  +  a12Ui 


xy 


>1 


+  (^4  —  x"bio)U2x  +  b8U2xx  +  bi2U2xy  +  Ai  +  0(h), 


T?j  =  Vec- 


^2  +  (C4  —  x"cw)Ulx  +  CsUixx  +  C\2Ui, 


xy 


(5.16) 

(5.17) 


+  (d^  —  x"  dio)U2x  +  d8U2xx  +  di2U2xy)  +  A2  +  0(h), 
where  Vi  and  V2  are  twelve  dimensional  vectors  similarly  defined  as  U±x,  Uyxx ■  U\xy, 
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U'2x-,  U‘2xx  ,  U-2xyi  and 

Vi(l)  =  ai  +  02, 

Vl(3)  =  05  +  06, 

foi  (5)  =  a9  +  aio  -  (sin2  0  +  1  - 

Vi(7)  =  b1  +  b2, 

Vi(9)  =  65  +  be, 

Vi  (11)  =  69  +  0io  +  sin0  cos  0, 
fo2(l)  =  ci  +  c2, 

V2(3)  =  C5  +  C6, 

V2(5)  =  eg  +  cio  +  sin  0  cos  0, 
V2(7)  =0i+  02, 

V2(9)  =  05  +  06, 

+2(11)  =  09  +  0io  “  (cos2  0  +  1 


Vi  (2)  =  a3  +  x"aio, 

Vi  (4)  =  07  —  (cos2  0  +  1  —  2  u~) 
2u~),  Vi  (6)  =  an  +  2  sin  0  cos  0, 

Vi  (8)  =  03  +  x"6io, 

Vi  (10)  =  07  —  sin  0  cos  0, 

Vi  (12)  =  0ii  —  (cos2  0  —  sin2  0), 
Vi  (2)  =  c3  +  x"cio, 

V2  (4)  =  C7  —  sin  0  cos  0, 

Vi (6)  =  cn  —  (cos2  0  —  sin2  0), 

Vi(8)  =  03  +  x"dvh 
Vi(10)  =  07  —  (sin2  0  +  1  —  2u~) 
-2v~),  Vi(12)  =  0ii  —  2sin0cos0, 


Ai  —  (04  —  aio  •  Xr,(0))«;i^o  +  a8wi^0  +  ai2'w1^Vo 

+(04  -  0io  •  x"(0))w2£o  +  b8w2^0  +  012^2^0  -  C\3, 
A 2  =  (c4  -  CIO  •  x"(0))wi£o  +  C8W  1^0  +  Cl2Wl^0 

+  (04  -  010  •  x"(0))w2^o  +  d8W2^0  +  012^2^0  “  Clj- 


where  Vl(i),  V2 (i)  are  the  i-th  components  of  El  and  V2  respectively. 

To  minimize  the  truncation  errors,  we  choose  the  coefficients  a%,  /3  k,  7 &,  and  +  so 
that  the  coefficients  of  ,  u^.  uUj,  u^.  uUjrj,  u^.  u^,  u~^,  u^,  u~^,  u^TjT] ,  and  u~^ 
vanish  in  (5.16)  and  (5.17).  Hence  we  set  the  following  system  of  equations  by  (5.16) 
and  (5.17): 


Vq  +  (04  —  x"aio)Uix  +  a8U\xx  +  a12Ulxy 

+  {b4  —  x"bio)U2x  +  b8U2xx  +  0i2 U2xy  =  0  , 

^  +  (C4- X^Clo)  fix  +  C8Ui  XX  +  C\2U\Xy 

+  (04  x  d\o)U2x  +  d8U2xx  +  d\2U2xy  =  0  , 
Ai  =  0;  A2  =  0. 


(5.20) 

(5.21) 

(5.22) 


The  system  of  equations  (5.20)  and  (5.21)  can  be  solved  separately.  Notice  that  in 
(5.20)  or  (5.21),  there  are  18  unknowns  and  12  equations.  The  optimization  method 
(cf.  [6,  17])  is  used  to  solve  the  system  of  equations.  Once  the  coefficients  a*,,  /3k,  7 k, 
and  8k  are  obtained,  the  correction  terms  C\ 3  and  Cf3  can  be  obtained  directly  from 
(5.22)  as  follows: 

C}j  =  (a4  -  a  10  •  x"(0))wi£o  +  a8w^0  +  a12wlivo 

+(b4  -  0io  •  x"(0))w2?o  +  08W2C&  +  bi2w2^0, 

Cfj  =  (C4  -  Cio  '  X"(0))W1^0  +  C8W^^0  +  Ci2Wizm 

+  (04  -  010  •  x"(0))w2£o  +  d8W2££0  +  d12w2^0. 

To  make  the  system  of  finite  difference  equations  better  conditioned,  we  impose  the 
sign  restriction  on  the  coefficients  ctk  in  (5.20)  and  the  coefficients  8k  in  (5.21) 


0:5  <  0  if  k  =  5, 


ctk  >  0  otherwise, 


(5.24) 
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c>5  <  0  if  k  =  5,  8k  >  0  otherwise. 


(5.25) 


We  also  impose  the  following  restrictions  for  all  the  coefficients: 

-1.1  *  bnd  <  ak,fik,lk^k  <  1-1  *  bnd  (k  =  1,2,  ••  •  ,9), 


(5.26) 


where 


bnd  =  max 


2(3  —  4u~)  2(3-4i/+)  1 


(5.27) 


With  these  restrictions,  the  solution  for  akl  /3k,  7/c,  and  8k  obtained  by  the  optimization 
method  are  of  order  0(l/h2). 

Once  we  have  all  the  coefficients  in  (2.1),  we  can  set  up  a  large  system  of  linear 
equations  (2.2),  (5.1)  and  (2.6),  (2.7)  with  2 (to  —  l)(n  —  1)  unknowns.  To  solve  such 
a  system  of  linear  equations,  we  have  tested  SOR,  preconditioned  GMRES(m),  and 
Bi-CGSTAB  methods.  Comparisons  among  these  methods  show  that  the  precondi¬ 
tioned  GMRES(m)  [4,  8]  and  preconditioned  Bi-CGSTAB  [4,  7,  14]  are  almost  equally 
successful.  We  will  give  some  results  in  the  next  section  using  the  GMRES(m)  method 
with  a  diagonal  preconditioning. 


6  Numerical  Examples 

We  have  performed  a  number  of  numerical  experiments  on  Sun’s  Ultra-10  workstations. 
In  these  numerical  experiments,  the  computational  domain  is  a  rectangular  region  with 
either  an  ellipse  or  a  circle  interface  within  the  domain. 

We  present  two  numerical  examples.  We  define  the  following  quantities 


Halloo  =  max  1 1  Ui(i,j)  -  Uuj\\ ,  Ratioi  = 

i-j 

||42)lloc  =  max  1 1 u2(i,j)  ~  U2ij\\,  Ratio2  = 


Halloo 

Halloo’ 

Halloo 


They  are  used  in  Table  6.1  and  Table  6.2. 

Example  1.  In  this  example,  the  domain  is  ST  =  (—1/2,12)  x  (—1/2, 1/2)  and  the 
interface  is  defined  by  x2  +  4 y2  =  0.352.  The  values  of  the  Poisson  ratio  and  the  shear 
modulus  are,  respectively, 


u+  =  0.20 
u~  =  0.24 


H+  =  1,500,000  in 

y~  =  2, 000, 000  in  ST 


The  Dirichlet  boundary  condition  and  the  interface  conditions  are  determined  from 
the  following  exact  solution: 


ui(x,y)  = 


U2  (x,y)  = 


xy  +  sin(l  +  x2  +  y2)  —  3  x2  +  y2 


in  fl+, 


xy  +  sin(l  +  x2  +  y2)  —  2x2  +  by2  —  4  in  G  , 

cos(l  +  x2  —  y2)  +  5  x2y  +  x2  —  y2  +  2  in 

cos(l  +  x2  —  y2)  +  5  x2y  +  3  x2  +  7 y2  —  6  in  ST 
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Table  6.1  gives  the  grid  refinement  analysis,  where  m  is  the  number  of  iterations 
before  restarting  in  the  GMRES  method.  Second  order  accuracy  is  achieved  since  the 
error  ratios  approach  number  four.  Figure  6.1  shows  plots  of  the  computed  solutions 
u\  and  U‘2  and  the  maximum  error  with  the  number  of  subintervals  along  each  side 
n  =  40  and  m  =  30. 


n 

1 1  -^n  \  \  oo 

Ratioi 

\\E^  II 

1 1  \  \  oo 

Ratio2 

m 

iterations 

20 

2.5749e-3 

* 

1.9612e-3 

* 

20 

174 

40 

5.0891e-4 

5.0595 

4.1027e-4 

4.7802 

20 

429 

80 

1.4667e-4 

3.4699 

9.4880e-5 

4.3241 

25 

1681 

160 

3.6265e-5 

4.0443 

2.3726e-5 

3.9990 

30 

2873 

320 

8.6890e-6 

4.1737 

5.4118e-6 

4.3841 

32 

8831 

Table  6.1:  The  grid  refinement  analysis  for  Example  1. 


GMRES(m)  method:  The  computed  solution  for  ul 
m  =  n  =40,  with  20  restarts 


GMRES(m)  method:  error  plot  for  ul 
m  =  n  =40,  with  20  restarts 


x  lO  4 


x  lO  4 


Figure  6.1:  The  computed  solution  and  the  maximum  error  of  Example  1. 
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Example  2.  In  this  example,  the  computational  domain  is  II  =  (—1, 1)  x  (—1, 1)  and 
the  interface  is  the  circle  x2  +  y2  =  1/4.  The  Dirichlet  boundary  condition  and  the 
interface  condition  are  obtained  by  the  following  exact  solution: 


ui  (x,y) 
u2{x1y) 


-(r4  +  c0  log(2r))/10  -  rg  +  (r$  +  c0  log(2r0))/10  in  12+, 
— r2  in  12“ , 

ln(l  +  x2  +  3 y2)  +  sin(a:y)  —  4r2  +  4cq  in  12+ , 
ln(l  +  x2  +  3y2)  +  sin(a:y)  in  12_, 


where  tq  —  0.5,  Co  =  —0.1,  and  r  =  \J x2  +  y2.  The  values  of  the  Poisson  ratio  and  the 
shear  modulus  are,  respectively, 

u+  =  0.20  in  fi+ 

v~  =  0.24  in  Q~ 

Table  6.2  gives  the  grid  refinement  analysis.  Second  order  accuracy  again  is  verified. 
Figure  6.2  shows  plots  of  the  computed  solutions  u\  and  U2  and  the  maximum  error 
with  the  number  of  subintervals  along  each  side  n  =  40  and  m  =  30. 


h  = 


H+  =  25,000,000  in  11+ 

=  30, 000,000  in 


GMRES(m)  method:  The  computed  solution  for  u  1 
m  =  n  =40,  with  25  restarts 


GMRES(m)  method:  error  plot  for  u  1 
m  =  n  =40,  with  25  restarts 


x  1  O'4 


GMRES(m)  method:  The  computed  solution  for  u2 
m  =  n  =40,  with  25  restarts 


GMRES(m)  method:  error  plot  for  u2 
m  =  n  =40,  with  25  restarts 


40 


Figure  6.2:  Computed  solutions  and  the  maximum  error  for  Example  2. 
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n 

1 1  -E'ra  1 1  oo 

Ratioi 

WE^W 

Ratio2 

m 

iterations 

20 

9.8376e-3 

* 

3.4529e-2 

* 

20 

158 

40 

1.9038e-3 

5.1675 

6.6401e-3 

5.2001 

25 

415 

80 

4.3297e-4 

4.3970 

1.5098e-3 

4.3979 

25 

1161 

160 

1.0680e-4 

4.0539 

3.6397e-4 

4.1483 

30 

3723 

320 

2.5572e-5 

4.1765 

8.5800e-5 

4.2420 

35 

8702 

Table  6.2:  Grid  refinement  analysis  for  Example  2. 


7  Conclusions 

In  this  paper,  we  have  developed  an  immersed  interface  method  for  elasticity  problem 
with  interfaces.  We  use  an  optimization  method  to  determine  the  coefficients  of  the 
finite  difference  equations.  We  employ  the  GMRES(m)  or  Bi-CGSTAB  to  solve  the 
large,  sparse,  and  non-symmetric  linear  system  equations  arising  from  the  discretiza¬ 
tion  of  the  elasticity  equation  together  with  the  interface  condition.  Our  numerical 
experiments  confirm  that  our  scheme  is  second-order  accurate. 
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